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NOMENCLATURE 


a(t) 

a 


c 

C. 


CN 

CY 



d 

D 

f 

z 


Two dimensional circular cylinder radius, — R(Z) 

2a 

Characteristic length of two dimensional unsteady flow, 
d/2 for constant afterbodies. 

Constant in Thwaites shape factor equation, = 5.1. 

sectional normal force coefficient. Normal Force Per Unit Length 

’ q d 

^00 

Sectional side force coefficient, 

q d 
^00 

Normal force coefficient, N/q^ S 
Side force coefficient, Y/q^ S 
Pitching moment about the nose, M/q^ Sz 

Pitching moment about (0, 0, \) 

Yawing moment about the nose> MY/q^ Sz 

Yawing moment about {0, 0, x) 

Three dimensional coefficient of pressure, P-P /q 

Two dimensional coefficient of pressure, p-p^/q 

Two dimensional drag coefficient, D/2?q 
Two dimensional lift coefficient, L/2aq 
Maximum body diameter 
Drag 

Fineness ratio, z/d 
Body length 


m 


M 

N 

P 

^co 

q 

R 


s 

t 

u 

V 

w 


[u. V] 

[u, V] 

[Ue , v^] 

[X, Y, Z] 


[x. y] 
[r, 0 ] 

GREEK 

ot 

r 

X 


Distance from cylinder surface at which vortices are 
introduced into the outer flow. 

Pitching moment 

Normal force 

Three dimensional pressure 

Free Stream dynamic pressure, l/2pV^ 

Cross flow dynamic pressure, l/2plj2 
Three dimensional body radius 
Vortex core radius 
Frontal area, iTd^/4 
Nondimensional time, Ut*/a^ 

Cross flow velocity, V sin a 
Fre& stream velocity 

Axial component of free stream velocity, V Cos a 
Boundary layer polar velocity components 
Cartesian velocity components 
Polar velocity components 

Polar velocity components at the cylinder surface 

Three dimensional coordinate system, Z positive along the 
axis from the missile nose. 

Two dimensional coordinate system 

Two dimensional polar coordinate system 


Angle of attack 
Circulation 

Moment anti, moment center distance from nose 

i 

• • • 

m 


Coalescence radius 


e 

p Density 

0 Vortex strength factor 

$ Potential function 

T Stream function 

0) Vorticity 

Sub- and Superscript 
( )* Dimensional 


Point vortex n 


1.0 INTRODUCTION AND SUMMARY 


The increased performance requirements of current aircraft and 
missile designs has caused a resurgence of interest in the aerodynamics 
of bodies at high angles of attack. In the past most aerodynamic 
design considerations were for angles of attack less than about 30 
degrees where either linear thoery applied (a < 5-7®) or the vortex 
formations due to boundary layer separation were symmetric and could 
be adequately handled using available empirical theories. New missile 
and aircraft designs, however, require consideration of angles of 
attack up to 180 degrees. Asymmetric vortex development for angles of 
attack greater than about 30 degrees require techniques for predicting 
the resulting side forces and moments. 

Once flow separation occurs on the lee side of an aerodynamic 
body the assumptions of inviscid flow and linearized theory are invalid 
and the only general approach is to solve the Navier-Stokes equations. 

The numerical solution of the three dimensional Navier-Stokes equations 
is not at hand for general missile or aircraft shapes. 

Other approaches to the problem must be sought for the present. 

One alternative to solving the full equations of motion entails the use 
of the viscous cross flow analogy. It was noted by early investigators 
that for angles of attack where the flow is assumed to be steady, if 
the flow is viewed in a cross flow plane, it appears similar to that of 
the time-development of flow behind a two-dimensional circular cylinder 
started impulsively from rest. This similarity is the basis of the 
viscous cross flow analogy in which the flow field and forces and moments 
on a missile are predicted from the analogous two dimensional time- 
dependent cylinder flow. Virtually all of the existing methods for 
predicting the aerodynamic characteristics of missile configurations at 
high angle of attack are based on the use of the viscous cross flow 
analogy. The methods have met with varying degrees of success. 

On the simplest level, the steady state value of cylinder drag 
suitably corrected for Mach number and Reynolds number effects is used 
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to compute the normal force distribution. This procedure was initially 
proposed by Allen (1) but has been extended by Jorgensen (2, 3, 4) and by 
Saffel (5) who have applied it from subsonic to supersonic flow for 
bodies of noncircular cross section, for angles of attack up to 180 degrees, 
and for wing-body-tail configurations. Thompson (6) has taken the 
analogy one step further and uses Sarpkays's (7) experimental impulsive 
flow circular cylinder data to represent the viscous contribution to 
normal force. Correction tables allow for Mach number, effect of nose 
shape, and transition from laminar to turbulent flow to be accounted for. 

The most recent cross flow models are those which replace the 
separated shear layer with a number of free point vortices which are 
allowed to roll up to form a concentrated vortex. Angelucci (8) and 
Marshall and Deffenbaugh (9) used this model for symmetric shedding. 

Wardlaw (10) and Deffenbaugh and Koerner (11) have extended the work 
to include asymmetric vortex development. 

If a large number of point vortices is used to model the wake 
computer time can become prohibitive. Only a limited number of test 
cases have been run using a large number of point vortices at high 
angles of attack for medium fineness ratio bodies (f 10.0). 

The objective of the current program is to modify the discrete 
vortex wake method of Reference (11) to efficiently compute the aerodynamic 
forces and moments on high fineness ratio bodies. The approach is to 
increase computational efficiency by structuring the program to take 
advantage of new computer vector software and by developing new 
algorithms when vector software can not efficiently be used. 

An efficient program was written and substantial time savings 
achieved. However, greater time savings were realized by using efficient 
FORTRAN programming techniques and by developing new algorithms than by 
vectorization. Several test cases were run for fineness ratios up to 
f = 16.0 and angles of attack up to 50 degrees. Prior to running the 
test cases a parametric investigation was carried out on an ogive nose 
fineness ratio, Ji^/d = 2.598, afterbody fineness ratio, i/d = 7.754, 
body. The parameters varied involved numerical and empirical parameters 


2 


used in the discrete vortex wake method. The results clearly point 
out limits of the method and its advantages and will be useful to 
anyone using this or similar methods involving discrete vortex wake 
approximations. 


2.0 TIME DEPENDENT VISCOUS CROSS FLOW 


For aerodynami c bodies at small angles of attack, the boundary 
layer remains attached except in the base region. For this situation, 
the boundary layer may be ignored allowing for an inviscid flow 
approximation and for slender bodies, application of linearized theory. 

For higher angles of attack regions of vorticity caused by flow separation 
along the body form on the lee side. Depending on the angle of attack 
these body vortices may be symmetric or asymmetric. Asymmetric vortex 
development can lead to large side forces and is currently the problem 
of most interest. 

An empirical method known as the cross flow analogy has been in 
use for several years and forms the basis for most of the high angle 
of attack prediction techniques used today. The difference in cross 
flow methods lie in the determination of the cylinder lift and drag 
coefficients. In the present study the time dependent flow about the 
circular cylinder is determined by coupling an inner boundary layer 
solution with an outer potential flow solution which models the wake 
using a large number of discrete vortices. 

2.1 CROSS FLOW TRANSFORMATIONS 

The geometry transformations for the cross flow analogy are obtained 
by considering the transformation 

(X*, Y*, Z*, T*) = (x*, y*, z* + W t*, t*} (2-1) 

and considering the flow in the z* = 0 plane. The basic cross flow 

assumption is that the flow in the X*, Y* plane of the body at some Z* is 

equivalent to the flow development at time t* around a two dimensional 

circular cylinder of radius where 

t* = Z*/W (2-2) 

(from (2-1) taking z* = 0.). 
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If the following nondimensional variables are defined 

Z = l*/i, a = a*/a, R = 2R*/d. t = ^* 

a 

where U is the cross flow velocity ( U= V cos a), i is body length, d 
is maximum body diameter, and a is a characteristic length of the two 
dimensional problem, then the basic cross flow equations become 


t = Z 4 tan a 

(2-3) 

a 

a(t)=|R(Z) 

(2-4) 


For most missile shapes a may be defined as d/2. However, for some 
shapes such as an ellipsoid which have decreasing afterbody diameters 
it may be more convenient to choose a as an average of the body radius, 

a = ^ J R* (Z*) dZ*. 

0 

2.2 FORCE TRANSFORMATIONS 

The normal force, N, and pitching moment, M, on the three dimensional 
body of revolution are to be obtained by integration of the normal force 
per unit length dN/dZ* over the body length, where dN/dZ* is assumed to be 
the same as the cross flow drag D(t*) on the two dimensional cylinder for 
t* = Z*/W. Similarily the side force, Y, and yawing moment, MY, are 
obtained by integration of the side force per unit length dY/dZ*, which 
is taken as equal to the cylinder lift, L(t* = Z*/W). 

The normal and side force acting on the body are then, 
i 

N =/ dZ* (2-5a) 

0 

i 

V =/ ^* dZ* ■ (2-5b) 

0 
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. and the pitch and yaw moments about some moment center (X*. Y*, Z*) * 

(0, 0, xsl) are: 

£ 

MY^ = r (Z* + R* dZ* + XtY (2-6b) 

If the coordinates are non-dimensional i zed as in Section 2,1 and the 
forces and moments are non-dimensional i zed as: 


Cd 

_ dN/dZ* 
l/2p V2d 

(2-7a) 

C, (Z) 

_ dY/dZ* 
l/2p V2d 

(2-7b) 

CN 

N 

l/2pV2s 

(2-8a) 

CV 

Y 

l/2pV2S 

(2-8b) 

*•111, 

M 

l/2pV2Si 

(2-9a) 


= c + xcn 
% 

(2-9b) 


MY 

l/2pV2S£ 

(2-9c) 

C 

n 

= C„ + ACY 

(2-9d) 
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where S is the frontal area S = 7rd^/4, then the force and moment 
coefficients become. 


CN 



Cd 


(2-lOa) 



C . dZ + XCN 
dZ ° 


“) C„ dZ + XCY 
dZ ^ 


(2-lOb) 


(2-lla) 


(2-llb) 


where f is the fineness ratio, i/d. 

The time dependent lift and drag on the two dimensional circular cylinder 
are: 



^27T 

p* cos 0 a*(t*) d 0 + I T* sin 0 a* 

*0 


(t) d 0 


(2-12a) 

-2tt ^27t 

L(t*) = -I p* sin 0 a* (t*) d 0 + I t* cos 0 a* (t) d 0 
‘'o ■'0 

(2-12b) 


Introducing the non-dimensional quantities, 

- n* p* 




P = 


1/2 p U2 


1/2 p U2 


IV..T = 


1/2 p U2 



and defining lift and drag coefficients as, 

D 


^D = 


1/2 p U2 2a 


(2-13a) 


^L = 


1/2 p U2 2a 


(2-13b) 


gives, 


2rr 


'D = 1/2 \ 0 a (t) d 0 + 1/2 J T sin 0 a(t) d 0 

0 ^ 

(2-14a) 


/2ir 27T. 

C, = - 1/2 I C sin e a(t) d 0 + 1/2 I x cos 0 a(t) d 

^ Jn 


(2-14b) 


Now, cross flow analogy assumes that at t* = z*/W 


dZ* 


= D (t*) 


(2-15a) 


— = L ft*l 
dZ* ^ ^ 


(2-15b) 


or in non-dimensional form the basic cross flow analogy force transforma- 
tions are 


(2-16a) 

(2-16b) 



If the three dimensional pressure coefficient is defined as 


S3D = - P: 


(2-17) 


1/2 p 
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and shear stress is neglected, then the sectional normal force coefficient 
is 


(Z) = 1/2 j Cp 3 p cos 0 R(z)d e (2-18) 

The two dimensional analog, neglecting shear and rewritten here for 
convenience is 

cos e a{t) d e (2-19) 

Now since the drag and sectional force coefficient are related as in 
equation 2-1 6a the correspondence between two dimensional unsteady and 
three dimensional pressure is 

SsD = f s <2-20) 

The two dimensional unsteady flow solution is given in the next 
section where it is shown that with the present formulation the above 
relationships between pressures differ by a additive constant when 
a(t) is varying, ie., on the missile nose. 




3.0 DISCRETE VORTEX WAKE (DVW) CROSS FLOW 


The problem is to find the flow field induced by a circular cylinder 
of time-varying radius in a uniform stream of a viscous incompressible 
fluid. The time dependent lift and drag on the cylinder will be applied 
as in Section 2.0 to predict the forces and moments on a three dimensional 
body of revolution. 

The flow field is that of a uniform flow with no body present for 
t < 0. At t = 0 a circular body appears at the origin and then grows in 
time, and in most cases (corresponding to typical missile shapes with 
cylindrical afterbodies) finally reaches a constant value. During the 
initial stages a boundary layer is formed on the cylinder. At later 
times the boundary layer separates from the cylinder feeding vorticity 
into the wake. Initially two symmetric regions of vorticity form 
behind the cylinder. Eventually these vortices become unstable, one 
breaks away and flows downstream initiating the Karman vortex shedding 
process. 

3.1 PROBLEM APPROACH 

The basic solution technique is to treat the problem of the 
impulsively started cylinder by assuming that for high Reynolds numbers 
the flow may be divided into a region of viscous inner flow near the 
cylinder described by a boundary layer and a rear shear layer, and a 
region of essentially inviscid outer flow elsewhere. The outer flow 
consists of the usual potential flow about a circular cylinder plus 
a wake region. The wake region is modeled by adding to the usual 
circular cylinder potential flow (vortex, uniform flow, and doublet) 
the potential flow induced by a set of point vortices. The point 
vortices represent the vorticity introduced into the wake from the 
separating boundary layer and rear shear layer. Since ideal point 
vortices have infinite velocities at their centers it is necessary to 
include in the model some kind of vortex core inside which the ideal 
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point vortex velocity is modified. This then is the discrete vortex 
wake (DVW) model . 

In the following sections the basic equations for the discrete 
vortex wake model are presented. Except for the unsteady integral 
momentum boundary layer solution, most of the theory presented has 
been previously detailed in Ref. (9). The equations are briefly 
rewritten here for completeness and to accompany the computer program 
users manual which is a separate volume of this report. 

3.2 OUTER FLOW 

The outer flow is defined as the inviscid irrotational flow 
outside of the boundary layer and rear shear layer. The exact equations 
which govern the outer flow are 


9(i) , 

8t 

o 

II 

3 

> 

• 

(3-la) 

U =1( 

X V¥ 

(3-lb) 

= 

w 

(3-lc) 


The discrete vortex method approximates a solution to equations (3-1) 
by replacing the continuous vorticity distribution w by a finite sum 
of N delta functions such that 


N 

0 )(x) = 2 6(x - x^(t)) (3-2) 

The velocity field is then found by summing the velocity field of the 
individual point vortices. To account for the presence of the 
circular cylinder, image vortices are included. Vortex motion is 
accomplished by integrating the equations 


n 


(3-4) 


dt 




where is the velocity induced by the distribution (3-2) at the 
discrete point vortex location x^. The self induced velocity of the 
point vortex is neglected. 

Of particular importance is the relationship between the solution 
of Equations (3-2, 3-3, and 3-4) and the exact solution of Equations (3-1). 
As pointed out by Saffman (12) it is not clear that the solution of the 
discrete vortex method converges to the exact solution for fixed N as t 
Although, as Saffman mentions, Dushane (13) proved under restrictive 
conditions that a form of the discrete vortex method converges in some form 
to the solution of Equations (3-1) for fixed t as N ->«>. 

For practical uses of the method, i.e., fixed N and t, confidence 
in the method depends on comparison with exact solutions and in some 
cases comparison with experiment. 

The complex potential for a circular cylinder with time dependent 
radius in the presence of a set of discrete vortices is 


2 IT 

w = $ + i4/ = 2+ ^- aS£nz + An (— ) 
z Z-n 'a 


N 

^ ^ 2 7r r 

n = 1 


'’n 


a^z. 


(3-5) 


z- 


where 


z = X + iy = re^® 

z - z^ = (x - x^) + i(y - y^) = r^^e’®ln 


(3-6a) 

(3-6b) 

(3-6c) 
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z - 


Z P 
n 


a^x 


1 

X - " 

+ i 

y - " 

0 2 

^n 


1 

<M 

C 

: 1 


= r2„e’®2n (3-6d) 


On the cylinder surface z = a, r^^ = ^r-j^/fi,^, and 'i' = 0. 

The coordinate system is shown in Figure 3-1. The cartesian components 
of velocity (u, v) and the tangential and radial velocity components 
(Ug, v^) may be obtained by differentiating (3-5) 


* = -u + iv = e-’9 (-v^ + lUg) (3-8) 

The cartesian components are 

aZ(x^ - y2 _ 2ixy) aa(x - iy) j. i^o (x - iy) 

-U + IV = 1 ; ^ 


5 r (y - yn) + - Xn) 

+ [ — I 

n=l 27 t j-2 

in 


27T 


2 1 
a'^y 


f 2 1 

a^x 

''n 


n 

y ~ 

+ 1 

X “ 

A2 


Jl2 

L nj 


L n 


'2n 




(3-9) 


The polar velocity components evaluated at the cylinder surface 


% ^ = a. 0) 


Vo = v^ (r = a, 0 ) 


become. 


-v„ + iu„ = 
0 0 


-a + i 2 sin 0 + 


ir 

c 

2TT 


+ i 






ar' 


H.1 

a 


} (3-10) 


In 
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FIGURE 3-1 COORDINATE SYSTEM | 
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For completeness equation 3-5 includes a circulation term r^. This 
is essentially an additional point vortex placed at the cylinder center. 
For bodies with sharp trailing edges is usually chosen to satisfy 
the Kutta condition. For the case of the circular cylinder where no 
sharp edge exists other physical considerations should determine the 

value of r . 

0 

However, to date there is no clear consensus on how to determine r 

0 

for this problem. The approach in the present work is to choose = 0.0, 
but to leave the capability in the computer program for future use. 

Pressure 


The surface pressure coefficient is determined from the time 
dependent Bernoulli equation. In non-dimensional form the Bernoulli 
equation is 

2$.t + Ug + v^ + p = 2F(t) (3-11) 


where from (3-5) 

$ = COS0 jr + ^ 


a a £n r 




such that 


N 


= 2aa cos e d(aa) ^ ^ - E !jl “ Qpn' ®nl (3-13) 

dt 27 t 2tt ' ' 


♦■t - r 


where 


'In 


In 


'2n 


t dt 




y-a^y 


dt K 


n I J_ 

f r2 


2 n 


0.. = 


^n"n 




(3-14) 
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If the pressure coefficient on the surface Cp = p - 


and from (3-11) 


p^ = 2F(t) - 1.0 - lim 2$,^ (r, 0, t) 

r -*■ <» 


(3-15) 


where 


$,^(r -> «■, 0, t) = - ^ (aa) lim Jin r ^ 0 ^ 


The pressure coefficient becomes 


- u^ + 1.0 + 2 ^ (aa) lim Jin r 

^ p -> 00 


(3-16) 


(3-17) 


For times where a 0, ie., on the missile nose, the last term in 

(3-17) precludes an absolute determination of Cp. The computations of 

C_ in the present program do not include the Tog terms in (3-17). On the 
P • 

nose then, C is only qualitative, on the missile after body (a = 0) C 

r r 

is computed exactly. 

Lift and Drag 

In the present program lift and drag are obtained by numerically 
(Simpson's Rule) integrating the pressure coefficient distribution around 
the cylinder. The additive constant in Cp for a 0 is integrated out 
and does not effect lift and drag. Lift and drag due to shear stress are small 
compared to the pressure forces and are ignored in the present formulation. 

A consequence of using Equation (3-17) to compute the pressure is 
the occurrence of spikes in the distribution due to discrete vortices 
being too near the surface point at which the pressure is being calculated. 

A core radius is employed to account for the singular nature of 
the point vortex approximation. In the present model point vortices 
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within a radius, r^ = .05, of the point at which the pressure is being 
computed are ignored in the computation. Even with this use of a 
core radius, spikes in the pressure coefficient distribution may 
occur. The discrete vortices, especially the ones introduced from 
the rear shear layer, may tend to cluster in regions of low velocity 
on the rear of the cylinder. This concentration of vorticity which 
is in effect a small secondary vortex will cause a spike. Also in 
examining the results, it appears that the time derivative of the 
potential term is a larger contributor to the pressure than is the 
surface velocity. 

Experimental pressure data is usually time-averaged and any 
sharp pressure peaks in the data are averaged out. 

Instantaneous pressure data. Ref. (14), may exhibit pressure peaks 
which as in numerical models. may be a consequence of small eddies in 
close proximity to the cylinder surface. The numerical results could 
be filtered or smoothed before the pressure is integrated. The 
difficulty which then arises is how much of the spike is "real" and 
how much is not. In the present work the computed pressure distri- 
bution is not smoothed. As a result the lift and drag may exhibit 
high frequency variations. See Figure 3-2. However, the forces will 
average out for frequencies on the order of the vortex shedding 
frequency. 

Vortex Motion 

Equation (3-4) is integrated using a Euler scheme to determine the 
development of the discrete vortex wake in time. Thus at each timestep, 
t|^, a new vortex distribution is calculated at tj^ ^ ^ = t|^ + Atj^ as 

+ l’ “ "("n- ■>'n’ 

V‘k + *k> • 
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2 TIME DEPENDENT LIFT AND DRAG OF 
A CONSTANT RADIUS CYLINDER 




On occasion when a discrete vortex is close to the cylinder 
surface the Euler scheme will allow the new discrete vortex position 
to be inside the cylinder surface. To avoid this each new discrete 
vortex position is checked and when found to be inside the cylinder 
is moved to a point (r^ = a + .001, (tj^)) just outside the cylinder. 

The integration step is then completed using only the tangential 
component of the original discrete vortex velocity. 


Vortex Birth 

Discrete point vortices are introduced into the flow at the 
boundary layer and rear shear separation points 9^, and respectively. 
The determination of the separation points is discussed in the next 
Section (3.3). The vorticity flux across the boundary layer at is 



(3-19) 


where the prime quantities denote inner boundary layer variables. 
The vorticity flux out of the boundary layer in a time-step, Atj^, is 
summed into a point vortex of strength 

It 


The point vortex is placed in the outer flow at the separation angle 
at some distance, m, from the cylinder surface. Images are created 
simultaneously to satisfy the condition of zero normal flow at the 
surface. In the present work, m is chosen so as to satisfy the no 
slip condition on the cylinder surface. The just born vortex and its 

images induce a velocity at the surface u. = With the requirement 

Trm 

that u. be canceled by the outer flow velocity then, 

J 


m 



At, u„ 
k 0 

2ir 


(3-21) 
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Point vortices are introduced from rear shear layer separation points 
in the same manner. 

The distance from the cylinder surface, m, depends on the timestep 
At. For large At, the point vortex may be placed outside of the boundary 
layer. In keeping with the boundary layer approximation it may be that m 
should be required to be within the boundary layer thickness. This 
restriction would limit the choice of At. The point vortex, however, 
could be required to remain within the boundary layer and the no slip 
condition relaxed. Because the point vortex is in the outer inviscid 
flow there is no real requirement for no slip to hold. 

Vortex Coalescence 

Even with the use of large scale computers and vectorized solution 
techniques, it may be desirable at some point to coalesce two or more 
vortices into one. In the present method, if the option is desired, two 
point vortices at (x- ,y. ), (x-s, y.) of strengths r., r. are coalesced into 

I I J J * J 

one discrete vortex of strength 

= r. + Tj (3-22) 


at position 



h “ yj ^ y.v’ 

(Ir^-l + iTjl) 


(3-22a) 


(3-23b) 
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when the two vortices are within some specified coalescence radius, e, 
ie., when 

(x,- - Xj)2 + (y. - yj)2 < e (3-24) 

Now, If 

(Xj. - x^.)/r and (y^. - y^.)/r 

where 

r = (x - x^.)^ + (y - y^)^ 
or 

r = (x - Xj.)^ + (y - yj)2 

are of order 0(e ), 

then it can be shown that the coalesced vortex yields a velocity field 

= (u^. + Uj.) (1+0 (e2)) 

Vc = (v^ + v^.) ( 1 + 0 (e2) 

As will be discussed later in the results section, coalescence of the 
vortices can make a significant difference in the early time development 
of lift on the cylinder even though the coalescence radius is no greater 
than e = .1 . 

Wake Vortex Stability 

The development of the wake will continue to be symmetric until some 
sort of perturbation is introduced into the numerical scheme so far 
outlined. The type of perturbation and magnitude of perturbation will 
have a large effect on the early flow development although an equivalent 
steady state will be achieved eventually. In the present work the method 
of perturbation is to reduce the boundary layer discrete vortex strength 
coming from one side of the cylinder for a short period of time 
early in the flow development. 
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Vortex Strength Parameter, a 


In the two dimensional theory all of the vorticity created at a 
solid body surface must eventually end up in the wake. Experimentally, 

Page (15) observed that only about 60% of the vorticity which is 
generated in the boundary layer is later measurable in the wake. It is 
generally felt that vorticity from the rear of the cylinder which is 
opposite in sign to that of the boundary layer coalesces with the 
boundary layer vorticity and accounts for Page's observation. Since 
rear vorticity is accounted for in the present model no further reduc- 
tion of vorticity should be required or perhaps permitted, at least 
for the two dimensional case. 

Por use in the cross flow analogy however a further reduction of 
vorticity may be necessary. In previous work. Ref. (9), good results 
were obtained by multiplying the boundary layer strength by a reduction 
factor a < 1.0 such that 

n j 

u^ 

where is the just born vortex of strength Fj = At ^ . 

The use of the factor a in the cross flow analogy is not yet 
fully understood. So far no theoretical justification for its use has 
been found. However, if one considers the problem of a three dimensional 
separation it is not reasonable to expect that all of the vorticity 
at the separation ends up in the cross flow plane. The vorticity vector 
in a three dimensional flow has three components. The justification 
for the use of the two dimensional vortex strength parameter a may be 
found in a better understanding of the three dimensional flow separation. 

3.3 BOUNDARY LAYER 
Governing Equations 

The boundary layer equations non-dimensionalized as in the outer flow 
and written in polar form for an expanding radius circular cylinder in 
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incompressible flow are 


m 3u . u au . V 3u 

a 


a(aU) aU 
at 


9(aU) + 
ax 


a^u 

ay2 


(3-25a) 


(o\ X M + X. = i 
^ ' q 2 8x ay a 


(3-25b) 


For a cylinder started impulsively from rest the initial conditions are 


(1) t = 0, u = V = 0 
and the boundary conditions are, 


(3-25c) 


(1) y = 0, u = V = 0 

(2) y u aU(x,t) (3-25d) 

Here x and y are coordinates along and normal to the wall. 

Equations (3-25) have in the previous work been solved directly using 

an implicit finite difference technique developed by M. G. Hall (16). 

The finite difference methods suffers from two aspects. The finite 
difference method requires a large amount of computer storage and 
computer time. Secondly, the finite difference method does not allow 
for oscillation of the separation angles. The discrete vortex wake 
solution requires only the determination of the separation angles at 
each timestep. Therefore in past work the finite difference method 
was replaced by a quasi -steady separation scheme developed by 
Stratford (17, 18). However, it was still necessary to use the 
finite difference solution for early times where the boundary layer is 
unsteady, and for times where a 0. 

The objective of the present work was to increase the efficiency of 
the discrete vortex wake method through efficient programming techniques, 
by using "vector" software where ever possible, and by developing new 
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algorithms. Therefore in the present work a method was developed to 
predict separation' on a changing radius cylinder based on the integral 
form of the unsteady momentum equations. This method based on the work 
of Thwaites (19, 20, 21) replaces the finite difference method. For a 
missile test case with ogive nose 2.59d and 7.727d afterbody at 45° 
angle of attack the unsteady finite difference solution executed on 
TRW's CDC6600 in 1207 CPU seconds. The new unsteady integral momentum 
solution executed in 363 CPU seconds. 

Unsteady Integral Momentum Formulation 

Multiplying (3-25a) by (aU - u), adding to (3-25b) and integrating 
with respect to y from o to ® yields the integral form of the momentum 
equations for time dependent radius. 


ii- + 

9t 




U2 

( 6 i + 262) + — 
a^ 


862 

~W 



au \ 
ay ^0 


(3-26) 


where U„ = aU 
0 


and 


61 = 




Displacement Thickness 


Momentum Thickness 


The last term on the left side of the equation (3-26) is analogous 
to a blowing term. Solutions to equation (3-26) based on the assumption 
of a Pohl hausen velocity profile, i.e., 
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(3-27) 


^ = an + bn^ + cn^ + 



may fail for strong blowing or suction because may be greater than 

1, which Pohlhausen profiles do not permit. In fact, for this problem 

a Pohlhausen solution was initially tried and did fail. 

To circumvent this difficulty a velocity distribution is chosen 

for which ^ must always be £ 1. Following the method developed by 
^0 

Thwaites (19), if we consider the velocity distribution given by 



or 


where conditions at the boundary are 

F(0) = 0 
F(l) = » 

and require that S be monotonic for 0 £ S < 1, then equation 3-28 
represents a velocity distribution which has no maximum of u/U^ other 
than unity. 

From the definition of momentum thickness it can be shown that 
the following condition must hold 

00 

f (2s - 1) F(s) ds = 1 (3-30) 


(3-28) 


(3-29) 
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F(s) can be chosen in any way such that the above equation (3-30) is 
satisfied. The momentum thickness is no longer a function of the shape 
of the velocity distribution since there are many choices for F(s)l This 
allows velocity distributions to be added simply and without imposing any 
condition on the momentum thickness, 62 , ie., (3-30) is satisfied. The 
idea behind Thwaites method is to add a Blasius profile and a suitable 
separation velocity distribution. The details are the same as presented 
by Thwaites (19) and the form of F(s) is 



where B(s) is the Blasius profile and x a parameter. From (3-28) and 
(3-32) the following derivatives at the wall are obtained 



(3-33a) 


(3f^) , 2cx3 
^0 ay^ 0 

where 

4 

j(x) = 4.5345 (1 - X^) [1-^(1- ^)] + CA (3-34) 

When A = 0, -^ = B(s), the Blasius profile. When A = 1, (|^) = 0 ie., 

02 ay 

separation occurs. For the above formulation to be complete a value for 
the constant, c, which appears in equations (3-32), (3-33), and (3-34) 
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needs to be determined. Thwaites (19) found that for the case of 
Howarth's linear adverse velocity gradient U = 3^ - 3ix, a choice of 
c = 5.1 resulted in separation of the boundary layer for 3i x/ 3 ^ = .120, 
ie*, the solution obtained by Howarth. 

Thwaites applied this method to a variety of steady problems 
Refs. (19, 20, 21) with suction. In the present work Thwaites method 
is extended to unsteady flows with a blowing like term. The value of 
c = 5.1 found by Thwaites is retained in the present work. 

Having chosen the form of the velocity profile as in (3-32) and 
defining the quantities 


H=5i/62, 2=6^. (|^)^ 

Then the integral momentum equation (3-26) can be written after a few 
pages of algebra as 



a^ 3x 




3Z 

3t 


4cx3 

- 2e + 


* 

( H + 2) + 2/T I 


-2z 


1 dH 


3( 1 - Xk) dx 
j 


X 3g 1 9U 

- — - 2 — — 

g at Uo at 


(3-35a) 


For the case of impulsive start the initial condition is, 

I.C. t = 0 ; 2 = 0 on X (3-35b) 

ie., momentum thickness is initially zero. 

In equation (3-35) g is a function only of the outer flow ie.. 
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9 = 

1 

a u 1 

0 + 

au 

0 



ax U 

0 

3t 

39 


1 1 

a^-u 

— 

= 

+ 


3t 


a^ 9 X. at U 

at^ 


.i- /!!^Y 

U \3tj 


(3-36a) 


2a aU (3-36b) 

0^ 

a^ 3x 


Now, from eqns (3-33) and evaluating (3-25a) at y = 0 


2cX^ 

±5^ = -2 g 


hence knowing z and g 
X can be found 

dH 

and H, e, j, and h are all functions of x, ie.. 


H = 2.5911 


^ _ c 
dx " 6 


1 - 


cx 


(1 -I-) 


, cx (x2 + 3) 
6 


e = 


62 


2.5911 (1 - X**) + x3(3X + 1) + 3 
1 - X^ 


{ ^ ) 

^ ay 0 


j(x) 


(3-37) 


(3-38) 

(3-39) 

(3-40) 


j(x) = 4.5345 ( 1 - X2) 
k(x) = 4.5345 { 1 - x2 ) 


1 - M - f ) 


( 1 - X- ) 


+ cx (3-41) 


- 2(4.5345)X 


1 - 


cx 


( 1 - ^- ) 


+ c 


(3-42) 
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differential equation, the 


ing the following set of ordinary 



(3-44) 

(3-45) 


(3-46) 


Equations (3-43) are integrated numerically using a simple Euler scheme. 
Once z is found as function of x and t then x is determined from (3-37). 
Separation occurs at x = when 

z(x = X , t) = — 

S nr2 


i e . , when A = 1 . 

Figs. 3-3 and 3-4 present the separation angles calculated using 
the method just outlined and using the previous finite difference 
scheme. The separation angle predicted by Thoman, (22) using a Navier- 
Stokes finite difference solution is presented for comparison. The 
results for a constant radius cylinder are presented in Fig. 3-3. The 
boundary layer finite difference and integral momentum methods are in 
close agreement. Figure 3-4 presents the results for a changing 
radius cylinder. (Thoman 's results are for constant radius). The 
results in Fig. 2-4 correspond to the rate of change of radius over 
an ogive nose of fineness ratio 2.598. A constant radius would be 
reached for this case at a time of t = 5.0. The time dependent radius 
results show that for early times separation is promoted over the 
constant radius case while for times t > 1.2 separation may even be 
prolonged. A possible explanation for this behavior might be found 
in examining the relationship between the -L terms and the blowing term 
- — U_ in equations (3-26) and (3-35) as, a, increases. The dependence 

o 0 

of the solution on the relationship between, a, and, a, has not been 
investigated in detail in the present work. 

Quasi -Steady Separation 

Stratford (17, 18) derived a method to determine the separation 
point of a laminar or turbulent boundary layer from an arbitrary pressure 
distribution. For laminar flow separation occurs at x when 
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Cp(x)| 


where the constant Sj^ = .14, and x is an equivalent distance, 

|.\ (U 

1 


X=X-X^-X +r'’’(77)^dx 

s ) u 


(3-47) 


which accounts for regions of positive pressure gradient. For 
turbulent flow separation occurs at x when 

(2Cp)5/4 = 1,06 3 (lO'^R)^''^^ 

dR 


(3-48) 


where R = Reynolds number based on local value of equivalent distance, 

X, and peak velocity, u^, or velocity at transition, u^^, whichever 
occurs later. The equivalent distance, x, accounts for an initial 
region of favorable pressure gradient and for a boundary layer which 

is initially laminar. The distance along the surface from e = 0 to the forward 
stagnation point is x^. 

X = X - Xp - Xj + 38.2(Up/u^^)’''® 


f (^=dx 

<^\f r\ 


5/8 


(x^^)3/8+ ^0 (^)3 (3_4g) 

Xtr 0 


where 


% = “m ' V ''m- Xn, " 


"o = V = “'V>’ '‘o = ^tr> ' V 


The empirical factor chosen by Stratford is 3 = .73 and is the value 
used in the present work. Mendenhall, (23), suggested that a factor 
6' = 3 sin a be used to account for 3D boundary layer effects. 
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4.0 VECTORIZATION 


Vectorization refers to the organization of individual steps in 
an algorithm into loop operations on arrays of data. Current vector 
processors exploit mainly pipelining rather than parallelism to 
achieve vector speedup over scalar operation. 

Pipeline processors consist of high speed arithmetic units which 
are segmented into a series’ of serial operations. The idea of vectori- 
zation is to perform operations on long vectors (such as C(I) = A(I) + B(I)) 
so that an element of the vector flows into the pipeline at each clock 
cycle. A clock cycle is defined as the maximum time it takes for an 
element to pass through any segment of the pipe. The idea is to keep 
all segments of the arithmetic unit continually busy on different 
elements of the vector. Once the pipe is full a result will be produced 
every clock cycle. In a scalar mode an element would have to traverse 
the entire pipeline to produce a result. In scalar mode then the 
operation C(I) = A(I) + B(I) would produce one result every 6 clock 
cycles (see Figure 4-1) whereas in vector mode 6 results would have 
been computed. 

The discrete vortex wake method involves operations on a large 
number of point vortices which are essentially vector operations. Most 
of the computation time is spent operating on these arrays. The 
present study was aimed at restructuring the current scalar code to 
take advantage of vector processing capability. 

Vector CDC7600 

When the CDC7600 is used with the FTN FORTRAN compiler it operates 
mostly as a sequential machine. However, by making efficient use of 
certain features of the machine hardware, ie., the functional units, 
and the instruction stack unit. the CDC7600 can exhibit a limited vector 
capability. A high speed software vector processing subroutine library 
was written by R & D Associates (24,25) for use on the NASA AMES CDC7600 
computer. When these subroutines are implemented on the CDC7600 it 
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FIGURE 4-1. ARITHMETIC PIPELINE OF TEXAS INSTRUMENTS 
ADVANCED SCIENTIFIC COMPUTER 
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displays characteristics similar to other vector processor machines, 
ie., TI's ASC and CDC STAR. 

The approach in the present work was to restructure the discrete 
vortex wake program into efficient FORTRAN code which then could be 
replaced by calls to the RDA vector routines. In the listing of the 
code the serial FORTRAN is left in as comment statements. 

Discrete Vortex Wake Vectorization 

Since conditional branches can not be handled in the vector method 
it was first necessary to divide the loops into two segments. One 
segment would contain all calculations required which were not dependent 
upon the conditional. In this segment a number of scalars were changed 
to vectors so that at the cost of a small amount of storage the vector 
routines could be used. 

In some cases the equations that followed the conditional could be 
expressed in a manner such that a vectorization could be used. Arrays 
containing either zeroes or ones were generated depending upon the 
conditional and were used in the calculations. Hence terms would be 
used or dropped depending upon the value in these arrays. In other cases 
this method could be used but it was so complicated that the vector code 
actually took longer to execute than the loop code vnth the conditional. 

At the beginning of the present program a constant radius two 
dimensional circular cylinder version of the discrete vortex wake program 
(VTX) was vectorized and timing runs made on the CDC7600. This version 
of the program used Stratfords criteria for separation. Timing runs on 
the program indicated that about 75% of the total time was spent in 
computing the motion of the vortices. The results of the timing runs 
are presented in Figure 4-2. For a time where there are about 800 vortices 
in the flow the optimized FORTRAN code resulted in an increase in speed 
of 35 percent. Vectorizing the optimized code resulted in another 14 percent 
increase. 

The cross flow discrete vortex wake program was then rewritten in 
structured FORTRAN code and then vectorized. Several algorithms ■ 
were rewritten and in the region of the missile nose a solution technique 
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FIGURE 4-2 VTX TIMING RESULTS 





based on the integral momentum form of the boundary layer equations was 
developed to replace the finite difference solution. A summary of the 
timing results is presented in Table 4-1. These results were obtained 
with the integral momentum boundary layer solution. 

The original cross flow program using the finite difference boundary 
layer solution was run for the NlBl missile test case at 45 degrees angle 
of attack. Using the RUNX compiler the program executed in 2694 CPU 
seconds on TRW's CDC 6600. The unoptimized program using the unsteady 
integral momentum solution executed in 2122 seconds using RUNX. The 
vecotorized integral momentum solution on AMES CDC 7600 took 190 seconds 
to execute using the FTN compiler. Allowing a factor of 5.0 for CDC 6600 
to CDC 7600 conversion and a factor of 2.0 to equivalence RUNX and FTN 
compilers the estimated time for the original finite difference program to 
run for this case on the CDC 7600 would be 269 seconds. The final vectorized 
cross flow program is then estimated to be approximately 30 percent faster than 
the original program. The cross flow program could not achieve the same 
increase in speed as did the constant radius test case (twice as fast as the 
original) because much less of the code could be vectorized. 

Results obtained using the vectorized Discrete VORtex Crossflow 
Evaluator (DIVORCE) program are presented in the next section. The user's 
manual for DIVORCE is a separate volume of this report. 
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TABLE 4-1 VISCOUS CROSS FLOW PROGRAM TIMING RESULTS 
NlBl (A^/d*2.598, f=10.325), a « 30° 


Cross Flow 

CDC6600 

CDC7600 

Program 

CPU SEC 

CPU SEC 

Unoptimized 

425 


Structured 

317 

60 

FORTRAN 



Vectorized 


53 
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5.0 METHOD APPLIED TO BODIES OF REVOLUTION 


The vectorized discrete vortex wake cross flow program was applied 
to several missile test cases of various fineness ratios with varying 
nose geometry and nose fineness ratios. Table 5-1 lists the test cases 
for which results are presented. Except for bodies NC (20° cone) and 
NPP (pointed paraboloid), the nose geometries were sharp tangent ogives. 

The results for cases NC, and NPP were non-dimensional ized using an average 
characteristic radius instead of the base diameter. The highest fineness 
ratio body tested was the ogive nose fineness ratio 5.0, afterbody fine- 
ness ratio 11.0, N3C configuration. 

Input Values 

In addition to body geometry several input values are required by the 
current method. These values are numerical parameters such as the 
integration timestep. At, and empirical constants like the vortex strength 
factor a. 

The numerical parameters are: 

At = integration timestep 
e = coalescence radius 

KRCOAL = timestep, t|^, at which to coalesce rear vortices 

KCOAL = timestep, tj^, at which to coalesce all vortices 

ZI = axial station at which solution is started ZI 0.0 

ZPERT = axial station at which a perturbation to newly born 
vortex strengths is added 

ZPEND = axial station at which a perturbation to newly born 
vortex strengths is removed 

The empirical constants are 

a = vortex strength factor 

Aa = vortex strength perturbation a' - a (KO + Aa) 

c = constant in Thwaites method = 5.1 (see section 3.3) 

3 = Stratfords turbulent separation criteria constant - .73 

Sk = Stratfords laminar separation criteria constant = .022 

Xtr = transition location on circular cylinder in radius 
measured from forward stagnation point .5236 
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Angle of Attack 

Nose Fineness 

Fineness 

Afterbody Body 

Reference 

Body 

a, DEG 

Ratio £ /d 
n 

Ratio f=£/d 

Diameter d( inches) 


NlBl 

10-50 

2.598 

10.325 

7.6 

26 

N4B2 

10-50 

1.633 

8.094 

7.6 

26 

N4B7 

45 

1.633 

7.0 

2.6 


N2C7 

45 

3.5 

7.0 

2.6 


N3C7 

45 

5.0 

7.0 

2.6 


NC 

10-50 

2.84 

2.84 

7.4 

27 

{20° Cone) 

NPP 

10-50 

2.84 

2.84 

6.0 

27 

(Pointed 

Paraboloid) 

N3C 

10-50 

5.0 

16.0 

2.6 

28 

N3C1 

35-50 

5.0 

12.0 

2.6 

28 

N2C1 

35-50 

3.5 

10.5 

2.6 

28 


TABLE 5-1 MISSILE TEST CASES 



The choice of these parameters and their effect on the solution is 
discussed in the following. 

5:1 NlBl RESULTS 

The predicted normal force coefficient on NlBl versus angle of 
attack is presented in Figure 5-1 and is compared with MX results. 

Reference (26). Past experience with the discrete vortex wake program 
has shown that good normal force results can be obtained with the vortex 
strength parameter, a = 0.6. However these results. Reference (9), were 
for angles of attack less than 30 degrees. For larger angles of attack 
only a very limited number of cases have been run. Figure 5-1 shows that 
good results are obtained up to 30 degrees for a = .6, but at higher 
angles of attack normal force is underpredicted. The bar on the last 
symbol denotes a range of values which may be obtained if during the 
solution the vortices are coalesced. This effect will be discussed in 
detail later. 

At 45 degrees angle of attack normal force coefficient was predicted 
using values of a = 1.0 and a = sin a = .707. The predicted values are 
higher since increasing vortex strength produces a stronger wake which 
increases pressure drag on the two dimensional circular cylinder. Figure 
5-2 shows that predicted pitching moment coefficients agree well with experiment 
to 40 degrees angle of attack with, a = .6, but drops off at higher angles. 

The distribution of normal force coefficient per unit length is presented 
in Figure 5-3. Increasing the vortex strength parameter results in a better 
prediction over the nose region but then overpredicts the load on the 
afterbody. 

Prior to exercising the discrete vortex wake program to predict 
side forces and moments on the goemetries listed in Table 5-1 a parametric 


42 



study was undertaken to determine the sensitivity of the method to the 
various numerical and empirical parameters listed in Section (5.0). In 
general normal forces and moments were not largely affected (except for 
o). However side forces and moments were nearly always affected. 

The parametric study was undertaken for the NlBl body since the 
largest amount of experimental data was available for that body. 

Effect of Vortex Strength Parameter a 

The effect of a on side force distribution is not well understood 
as it is for normal force. While eventually (in the steady state) 
stronger vortices would probably result in larger amplitudes of 
oscillating lift it is not clear how stronger vortices affect the initial 
vortex development. The side force coefficient on NlBl is presented 
in Fig. 5-4. The magnitude of the predicted coefficients bounds the 
values obtained experimentally. However, at the same angle of attack the 
side force coefficient may be positive or negative depending on the 
choice of a. Figure 5-5 shows the side force coefficient per unit length 
along the missile axis for three values of a. The development of the 
asymmetry is similar up to about half (Z = .5) the body length. In this 
example a = .707, produced the largest side force. The corresponding 
yawing moment coefficient is presented in Figure 5-6. 

Effect of Coalescence on Side Force Distribution 

To save computer time it may be desirable to coalesce two vortices 
into one if they are in close proximity. Virtually all of the discrete 
vortex methods have some type of coalescence scheme. Normal force is 
affected little by coalescence. Side force on the other hand is signifi- 
cantly affected as is shown in Figure 5-7. Vortex plots at Z = .804 
are presented in Figure 5-8 for the case of no coalescence and coalescence 
using the largest core radius of e = .25. From Figure 5-8 very little 
difference in the wake structure can be observed (aside from fewer vortices). 
In examining the pressure distributions at Z = .7516, Figure 5-9, however* 
it appears that small differences in the pressure occur and result in 


the significant differences in side force when integrated. 

The initial vortex growth is very sensitive to small perturbations. 
Coalescing the vortices introduces an error of only order at the 
time of the coalescence. However this small perturbation significantly 
affects the initial vortex growth and eventual asymmetry. At later times 
a Karman vortex street will always be established and the initial pertur- 
bation is not important. However for times pertinent to the cross flow 
analogy coalescence makes a difference as indicated in Figure 5-7 for 
a = 30^ and Figure 5-10 for a = 45°. 

For the same reasons described above the frequency of coalescence, 
or the time at which the first vortices are coalesced can alter the 
distribution of side force, see Figure 5-11. 

Effect of At on Side Force Distribution 

In the present work At is usually allowed to vary from At = .05 
near the nose tip where the cylinder radius is small and a small time- 
step is required to a constant value usually At = .125 at the nose 
afterbody junction. Figure 5-12 shows the effect of constant At = 

.125, .250. Letting At be variable had little effect. However, At 
= .250 significantly reduced the amplitude of the sectional side force 
coefficient. For 30 degrees angle of attack At = .0625 was also 
tried with results significantly different from the At = .125, and 
At = .250 cases. However the maximum amplitude of the sectional side 
force coefficient for 30 degrees was only 0.1. 

In addition to being the integration timestep. At also implicitly 
determines the distance m (see Eqn 3-21) at which vortices from the boundary 
layer and rear shear layer are introduced into the outer flow. Further 
work needs to be done in examining the effect of At. The discrete 
vortex wake solution may not be unique as indicated by many of the 
present results. For many of the cases run in this and other studies 
At = .125 has provided reasonable answers. For the rest of the results 
presented in this work At = .125. 
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5.2 OTHER TEST CASES 


Based on the study of the NlBl body, discrete vortex wake parameters 
were chosen as follows to investigate the effect of angle of attack, 
nose fineness ratio, nose geometry, and overall fineness ratio: 
a = .6 (unless otherwise specified) 

Aa = . 1 
11 = .01 

AT = .125 (variable .05 to .125 on the nose) 

ZPERT = .05 
ZPEND = .25 

KRCOAL - 30| vortices coalesced, ie., normally 
KCOAL = 150) a = 10, 20, 30 does not require coalescence 

N4B2 

The forces and moments on N4B2 (^^/d = 1.633) are presented in Figs. 
5-13 through 5-17. The normal force coefficient agrees well with the 
experimental data up to 45 degrees angle of attack. The pitching 
moment coefficient is somewhat higher than the experimental values for 
angles of attack greater than about a = 30 degrees. Examination of the 
normal force distribution shows that the predicted normal force is 
higher over the nose region than that measured experimentally. This is 
in contrast to the NlBl results which underpredicted nose distribution 
of local normal force coefficient. Predicted and experimental side forces 
are small except at 45 degrees angle of attack. The predicted results 
show a change in side force sign only at a = 40 degrees. Similarly the 
predicted yawing moment changes sign whereas the force data does not. 

The predicted force and moment results are of the same order of 
magnitude as the experimental values. A plot of |CY| versus angle of 
attack would show better agreement. However, the purpose of the present 
work is to examine and disclose the applicability of discrete vortex 
methods for the cross flow problem and therefore signed values are 
presented. 


|(unless otherwise specified) 
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N2C1 , N3C1 


The present method was applied to the data taken by Jorgensen 
Ref. (28) on ogive cylinders of nose fineness ratio = 3.5, and 
A^/d = 5.0, and afterbody ratio 7.0. A calculation was also done 
with afterbody fineness ratio 11.0. The purpose was to determine if 
effects of nose fineness can be predicted using the discrete vortex 
method. The normal force coefficient data in Figure 5-17 shows that 
the predictions are significantly lower than the experiment. Part of 
the discrepency is a Mach number effect. Figure 5-18 shows the 
predicted normal force coefficient with a = 1.0 to be in much better 
agreement. Predicted aerodynamic normal force center also agrees 
well with Jorgensens results. Figure 5-19. Jorgensen found maximum 
side force to occur at 35 to 40 degrees on N3C1 , Figure 5-20. Maximum 
predicted side force coefficient occurs at 50 degrees angle of attack 
and is only about one half of the maximum value obtained by Jorgensen. 
Yawing moment is presented in Figure 5-21. Side force coefficient 
results for N2C1 and for N3C (f = 16.0) are presented in Figure 5-22 
in addition to the N3C1 results. The N3C1 configuration (£^/d=5.0, 
f=12.0) resulted in the largest predicted side force. 

Effect of Nose Fineness Ratio on Side Force Coefficient Distribution 

To examine in more detail the effect of nose fineness ratio on 
asymmetric vortex develpment and side force distribution three bodies 
of overall fineness ratio f = 7.0, and nose fineness ratios, A^/d=1.666, 
3.5, and 5.0 were run at 45 degrees angle of attack. Discrete vortex 
wake parameters were ZI = .01, ZPERT = .05, ZPEND = .25, a = .6, 

A 0 = .1, At = .125, and were the same for the three runs. The vortices 
were not coalesced. The sectional side force coefficients Figs. 5-23, 
all are initially perturbed to positive values and then reach a maximum 
negative amplitude at about, Z = .4. The longest nose attained the 
largest sectional side force amplitude. The effect of nose fineness 
ratio on separation angle is presented in Figure 5-24. Separation is 
in general further windward for the highest fineness ratio body. The 
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short ogive tends to prolong separation. The jump in the separation curve 
at Z = .3 for caseN4B7is a result of changing from the unsteady boundary 
layer solution to Stratford's quasi -steady scheme. The separation angles 
for all three bodies approach = 85 degrees at about 70% of the body 
length. Figure 5-25 shows the vortex growth for the short, N4B7, and 
long N3C7 ogives. The N2C7 development was very similar to N3C7. The 
radii are to the same scale. Some differences are observable in the 
vortex wake. However, overall the vortex development looks very similar. 

While nose fineness ratio has an effect on the predicted asymmetric vortex 
develppment and on the resulting side force distribution it is not clear 
that these effects are any greater than the effect of changing one of the 
other parameters in the method such as. At, or coalescence radius or 
frequency, and probably has much less of an effect than changing a, or 
ZPERT and ZPEND. Figure 5-26 shows the effect of applying the vortex strength 
perturbation over different lengths of a fineness ratio A^/d = 5.0 ogive 
nose. The length scale in Figure 5-26 is non-dimensional i zed by the 
longest body length, = 16.378 inches. 

The effect of forebody geometry was briefly examined for the case of 
a sharp 20° cone and a pointed paraboloid (Ji^/d = 3.5, f = 3.5). The 
side force coefficient is presented for comparison with the available 
data. Forebody geometry had little effect on the predicted side force 
coefficient. 

Pressure Distribution 

One of the purposes of a complex method such as the discrete vortex 
wake program is to be able to predict surface pressures and the 
distributed aerodynamic load. The sectional normal force coefficient and 
sectional side force distributions on NlBl at 45 degrees angle of attack 
is presented in Fig. 5-28, and 5-29. 

The corresponding vortex development and surface pressures are 
presented in Figures 5-30 and 5-31. The numerical perturbation caused 
an asymmetric vortex development which resulted in a side force distri- 
bution opposite in sign to that obtained experimentally. Applying the 
perturbation to the opposite point would change the sign of 
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the predicted results. Rather than rerun the case the results can be 
viewed as in Figure 5-30. The usual nomenclature is y positive to the 
right. Figure 5-30 shows the vortex development on NlBl at 45 degrees 
angle of attack. Corresponding pressures are presented in Figure 5-31. 

At Z = .124 the experimental and predicted pressures are almost 
symmetric. 

At Z = .222 the experimental pressures are still symmetric 
whereas slightly asymmetric pressures are predicted at Z = .232. 

Pressure recovery on the rear of the cylinder is considerably lower 
for the experiment. At Z = .337 the experimental pressure begin to 
show an asymmetry. The theoretical pressures at Z = .352 are asymmetric 
but with lower pressures at 80 degrees as opposed to the experiment 
which exhibits the lowest pressure on the opposite side at 285 degrees. 

At Z = .352 the vortex plot Figure 5-30 shows the left vortex 
beginning to dominate. Hence, the lower pressures on the left half 
of the body { 0 < 9 < 180 degrees). At Z = .594, Figure 5-30, the left 
vortex has entrained fluid from the right side and has weakened. This 
allows the right vortex to increase in strength as the weakening left 
vortex begins to move downstream. The corresponding pressure Figure 
5-31 ( Z=.594 ) shows a large asymmetry with the lowest pressure on 
the right side where the stronger right vortex exists. The agreement with 
the experimental pressures, Z = .591, is very good. At the last station 
the left vortex is again dominant as shown in Figure 5-30, Z = .835. The 
predicted pressures at this station show the presence of a strong left 
vortex. The experimental pressure data. Figure 5-31 Z = .845 reflects 
a similar asymmetry. 
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NORMAL FORCE COEFFICIENT CN 



FIGURE 5-1. NORMAL FORCE COEFFICIENT ON NlBl (s, /d = 2.598, 
f = 10.352) vs. Angle of Attack 
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PITCHING MOMENT COEFFICIENT, C (x - .647) 


• MX, Run .309, Re^ = 1.27 x 10^, M » .4 
X MX, Run 312, Re, - 3.17 x 10^, M » .4 



FIGURE 5-2 NlBl PITCHING MOMENT COEFFICIENT WITH ANGLE OF ATTACK 
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LOCAL NORMAL FORCE COEFFICIENT 



FIGURE 5-3 EFFECT OF VORTEX STRENGTH PARAMETER a, ON NlBl, a = 45° NORMAL FORCE DISTRIBUTION 
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MX, Run 312, Re^ = 1.27 x 10®. M*.4 
MX, Run 309, Re^ == 3.17 x 10^, M-.4 

DIV0RCE(o * .6, at = .125) 

DIV0RCE(o » 1.0, at » .125) 
DIV0RCE(ff - .8, at - .125) 

DIV0RCE(o *.707, at - ,125) 




SECTION/ 



DIVORCE {at-. 125) 

• a-1.0, CN-6.329, CY-.247 

• 0-.8, CN-5.281, CY--.526 

□ a-. 707, CN«4.944,.CY--1.358 




FIGURE 5-5 EFFECT OF VORTEX STRENGTH PARAMETER, a, ON NlBl SIDE FORCE DISTRIBUTION, a 


YAWING MOMENT COEFFICIENT C„ (A = .647) 


A MX, Run 312, Re^ - 1.27 x 10®, M - .4 
□ MX, Run 309, Re^ - 3.17 x lO^, M * .4 

• DIVORCE(a ** .6, AT - .125) 

• DIVORCE(or « 1.0, AT - .125) 

• DIVORCE(a = .8, AT » .125) 

• DIVORCE(a = .707, AT * .125) 
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SECTIpNAL SIDE FORCE COEFFICIENT, C^ 



FIGURE 5-7 EFFECT OF COALESCENCE RADIUS ON NlBl SIDE FORCE DISTRIBUTION, a=30° 
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FIGURE 5-8 EFFECT OF COALESCENCE ON DISCRETE 
VORTEX WAKE DEVELOPMENT. Z=,804 
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PRESENT THEORY (aT « .125. o • .6) 
• NO COALESCENCE < 



FIGURE 5-9 EFFECT OF COALESCENCE RADIUS, e, ON NlBl PRESSURE DISTRIBUTION, a=30°, Z =.7516 




DIVORCE (at-. 125, a-. 6) 



FIGURE 5-11 EFFECT OF COALESCENCE FREQUENCY ON NlBl SIDE FORCE DISTRIBUTION, a=30“ 


SECTIONAL SIDE FORCE COEFFICIENT 





NORMAL FORCE CN 


A MX. Run 232, Red«1.27 x 10^, M«.4, PRESSURE INTEGRATION 
□ MX. Run 174, Red-1.27 x 10‘. M-.4. GRIT RING FORCE DATA 
O DIVORCE, (a-. 6. AT-. 125) 


50 a DEG 60 


FIGURE 5-13 NORMAL FORCE COEFFICIENT VARIATION WITH ANGLE 
OF ATTACK N4B2 (ii„/d = 1.666, f = 8.094) 





PITCHING MOMENT COEFFICIENT C ( X = .548) 


F 


X MX, Run 174. N4B2G. Re^ - 1.74 X 10*. H - .4 
• DIVORCE (o ■ .6, 4T - .125) 



10 20 30 40 50 


FIGURE 5-14 N4B2 PITCHING MOMENT COEFFICIENT VARIATION 
WITH ANGLE OF ATTACK 
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SIDE FORCE COEFFICIENT CY 


A MX, Run 232, Re^ » 1 .27 x 10«, M-.4, PRESSURE INTEGRATION 


O MX, Run 174, Re^ - 1.27 x 10*, M«.4,. GRIT RING FORCE DATA 
O DTV0 RCE(o * .6, M = .125) 



Figure 5-15 SIDE FORCE COEFFICIENT VARIATION WITH 

ANGLE OF ATTACK (N4B2, ii^/d=1.666, f=8.094) 
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YAWING MOMENT COEFFICIENT C (X = .548) 


X MX. Run 174, N4 B2G, Re^ - 1.74 x 10®. M - .4 
• DIVORCE(o « .6, AT • .125) 



FIGURE 5-16 N4B2 YAWING MOMENT COEFFICIENT 
VARIATION WITH ANGLE OF ATTACK 
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NORMAL FORCE COEFFICIENT CN 


□ JORGENSEN, REF.28, N2C1 {/^/d=3.5, f=10.5) 

0 JORGENSEN, REF.28, N3C1 (j^^/d=5.0, f=12.0) 

^ DIV0RCE(a=.6, aT=. 125), N2C (l^/d=3.5, f*10.5) 

A DIVORCE(<j=.6, aT=. 125), N3C (l^d-S.O, f=12.0) 
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FIGURE 5-18 EFFECT OF VORTEX STRENGTH PARAMETER, 
a, ON N3C1 NORMAL FORCE 
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□ O 


AERODYNAMIC FORCE CENTER, XACN/D 


RE 

O 2.2 X 10^ Jorgensen, Ref 28 . N3C1 
□ 4.3 X 105 ^ . 5 0, t - 12.0) M - .6 
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• DIV0RGE(o - .6, at - .125) 

10 r ^ DIV0RCE(a - 1.0, AT » .125) 



FIGURE 5-19 EFFECT OF VORTEX STRENGTH PARAMETER, a. 


ON N3C1 AERODYNAMIC FORCE CENTER . 
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FIGURE 5-20 EFFECT ON VORTEX STRENGTH PARAMETER, a, ON N3C1 SIDE FORCE 
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YAWING MOMENT COEFFICIENT CYN 




SIDE FORCE COEFFICIENT CY 


O JORGENSEN (/^/d«3.5, f-10.5), REF 28 
□ JORGENSEN (/.^/cl-5.0, f-12), REF 28 
A KEENER & CHAPMAN (/^/d-5.0, f*16), REF 27 
<$> DIV0RCE.,£^/«1'3.5, f-10.5 
S DIVORCE, ^yd«5.0. f-12.0 
4 DIVORCE,/ /d-5.0, f-16.0 
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FIGURE 5-22 EFFECT ON FINENESS RATIO 
ON SIDE FORCE COEFFICIENT 
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SECTIONAL SIDE FORCE COEFFICIENT, C^ 
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FIGURE 5-25 




FIGURE 5-25 CONTINUED 
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DIV0RCE(o-.6. Ao-.lo, AT-.125) 

CASE 350 (l^/d-5.0, 1-16.378 In, f»16.0) Zi=.05, Za=.25, Zo=.01, CY=.816, CN-6.904 
CASE 2004 {l^/d»5.0, 1-14.843 1n, f»12.0) Zj=.022, Z2=.109, Zo=. 00435, CY=-1 .803, CN=5.085 



SIDE FORCE CY 
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FIGURE 5-27 EFFECT OF FOREBODY GEOMETRY ON SIDE FORCE 
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SECTIONAL NORMAL FORCE COEFFICIENT, 





, SECTIONAL SIDE FORCE COEFFICIENT 



FIGURE 5-29 SIDE FORCE DISTRIBUTION, a = 45° 
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N1B1(Z=.594) 
FIGURE 5-30 CONTINUED 
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FIGURE 5-30 CONTINUED 
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FIGURE 5-31 CONTINUED 







6.0 DISCUSSION 


Effective numerical algorithms were developed to apply the discrete 
vortex wake cross flow method of Ref. (9) to high fineness ratio bodies 
at high angles of attack of 30 to 50 degrees. Because of computer time 
and cost limitations the method had been used previously only on a very 
limited basis to demonstrate capability for predicting asymmetric wake 
development. 

The required computational ability was achieved by (1) developing 
a method to solve the unsteady integral momentum boundary layer equations 
on a changing radius cylinder to replace the current finite difference 
solution, (2) restructuring the program to take advantage of Vector 
Processor Architecture in the vortex array computations, (3) optimizing 
the FORTRAN code which could not be vectorized, and (4) using a bigger 
computer. 

After converting CDC6600 (RUNX) run times to CDC7600 run times 
using the optimizing FTN compiler, the increase in execution speed 
of the vectorized program over the original program was found to be 
about 30 percent. A 30 percent difference in execution speed alone does 
not justify using the method now, whereas it was too costly before^ 

However, coupled with the increased speed of the CDC7600, over TRW's 
CDC6600, the method was able to be tested on various missile geometries, 
and the effect of numerical and empirical parameters examined. 

Only about 12 percent decrease in computation time was achieved 
using the vector processing instruction set. Larger time savings were 
not achieved primarily because of two reasons. The first reason is 
linked to the fact that discrete vortex wake methods must include a 
vortex core radius to modify the singular behavior of point vortices. The 
vortex computations involve different operations depending on whether a 
calculation point is within a core radius of any point vortex. This core 
radius check interrupts the flow of data into the pipeline and results 
in minimal vector performance. 

Probably 60% of the computation time is spent in the point vortex 
computations. The other 40% of the code could not be efficiently vectorized. 


90 



The FTN compiler optimizes CDC 7600 code very well. To gain speed over 
the very efficient serial code, a significant portion of the code may 
have to be vectorized. 

The results obtained using the vectorized code, and the serial code 
were found to be slightly different. To examine this further, the serial 
program was stopped with about 1200 vortices in the flow. The vector 
program and the serial program were then started using the same 1200 vortex 
positions and strengths. The results were identical out to about the 12th 
significant digit. However, if the runs were continued for some time, 
slight differences would develop. 

Predicted normal forces and moments in general agreed well with 
experiment. Increasing the vortex strength parameter, a, resulted in 
better agreement with Jorgensen's (28) M = .6 data. 

Probably the most significant aspect of the present work was the 
discovery of the extreme sensitivity of the side force development on the 
numerical and empirical parameters used in the method. Perhaps this 
result should not be surprising since experimental side force data is also 
extremely sensitive to perturbations in the flow. 

Sarpkaya (29), in a recent experiment using a vertical water tunnel 
measured the time dependent drag and lift force acting on impulsively 
started circular cylinder. Experimentally Sarpkaya found that in impulsive 
flow the development of lift is strongly dependent on the initial disturbances 
which give rise to asymmetry and subsequent vortex shedding. In some runs 
Sarpkaya found that the vortices continue to grow symmetrically and result 
in small lift forces while for repeated runs at the same Reynolds number 
the asymmetry develops rapidly and the lift force not only starts at smaller 
values of t but also reaches larger amplitudes. 

The predicted results of the present study indicated the same type of 
randomness in the early development of lift that Sarpkaya observed 
experimentally. 

It is possible that the discrete vortex wake problem is not well posed, 
ie., a small change in the initial conditions result in a large change in 
the answer. However, it may be possible to bound the problem. By varying 
the parameters in the method within reasonable limits, it may not be 
possible to obtain a value larger than a certain maximum value. 
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7.0 CONCLUSION AND RECOMMENDATIONS 


The viscous cross flow analogy was employed using a two dimensional 
discrete vortex wake cross flow s^olution to predict the flow field 
around high fineness ratio bodies of revolution at high angles of 
attack up to fifty degrees. A method was developed to solve the unsteady 
integral momentum boundary layer equations for the expanding radius 
circular cylinder started from rest in laminar flow. Since the integral 
momentum equation is valid for either laminar or turbulent flow, an 
extension of the present solution technique using a modification of 
Thwaites shape factor could be developed to extend the method to high 
Reynolds number turbulent flows. 

Predicted normal forces and moments agreed well with experimental 
data. Predicted side forces and moments fell within the range of the 
experimental data but were found to be extremely sensitive to small 
perturbations. The maximum value of predicted side force was generally 
about half of the maximum experimental value. The predicted maximum 
values also occurred at different angles of attack. 

The distribution of side force is mostly affected by the numerical 
perturbation applied to the method to induce asymmetric vortex develop- 
ment and by the vortex strength parameter. Increasing the vortex 
strength parameter always results in an increase in normal force coef- 
ficient, but not necessarily in side force coefficient. A physical 
explanation of the vortex strength parameter may be found by considering 
that the probability that only part of the total vorticity shed from a 
three dimensional separation line ends up in the cross flow plane. 
Experimental and theoretical investigations of the three dimensional 
flow separation on a typical missile shape nneds to be undertaken to 
determine the cross flow component of vorticity. 

Forcing a small asymmetric change in the strength of the newly 
introduced vortices for a portion of the missile nose has a significant 
effect on the side force development. Forcing the separation angles to 
be slightly asymmetric would be a similar perturbation mechanism. The 
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initial disturbance may be due to asymmetries which develop first in the 
wake (ie., free stream turbulence) and then feed back to the boundary 
layer. The physics of the initial instability needs to be better under- 
stood in order to determine how to model the initial disturbance. 
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